Mon. Not. R. Astron. Soc. 000, 000-000 (0000) Printed 2 December 2011 



(MN WF$i style file v2.2) 



New Light-Travel Time Models and Orbital Stability 
Study of the Proposed Planetary System HU Aquarii 

^ ;T. C. Hinse 1 ' 2 , J. W. Lee 1 , K. Gozdziewski 3 , N. Haghighipour 4 
O | C.-U. Lee 1 , E. M. Scullion 5 

' ^ , 1 Korea Astronomy and Space Science Institute, 776 Daedukdae-ro, Yuseong-gu, 305-348, Daejeon, Republic of Korea. (tchinselSgmail.com) 
Q • 2 Armagh Observatory, College Hill, BT61 9DG, Armagh, UK. 

' 3 Nicolaus Copernicus University, Torun Centre for Astronomy, PL-87- 100 Torun, Poland. 

4 Institute for Astronomy & NASA Astrobiology Institute, University of Hawaii, 96822 HI, USA. 

5 Institute of Theoretical Astrophysics, University of Oslo, 0371 Oslo, Norway. 



PL, 

6 



> 

\o 
o 
o 



- 1—1 

X 



2 December 2011 



ABSTRACT 

In this work we propose a new orbital architecture for the two proposed circumbinary 
planets around the polar eclipsing binary HU Aquarii. We base the new two-planet, 
light-travel time model on the result of a Monte Carlo simulation driving a least- 
squares Levenberg-Marquardt minimisation algorithm on the observed eclipse egress 
times. Our best-fitting model with \r = 1-43 resulted in high final eccentricities for the 
two companions leading to an unstable orbital configuration. From a large ensemble of 
initial guesses we examined the distribution of final eccentricities and semi-major axes 
for different \r parameter intervals and encountered qualitatively a second population 
of best-fitting parameters. The main characteristic of this population is described by 
low-eccentric orbits favouring long-term orbital stability of the system. We present 
our best-fitting model candidate for the proposed two-planet system and demonstrate 
orbital stability over one million years using numerical integrations. 

Key words: binaries: close - binaries: eclipsing - stars: individual: HU Aquarii - 
planetary systems: circumbinary planets 



INTRODUCTION 



5-H . In the p ast two decades, astrophy sical timing measurements has been used to infer the existence of multiple low- mass planetary 



objects. Wolszczan &: Fraill fl992) announced the first detection of a planetary system around the pulsar PSR1257+12. A two- 
planet sys tem orbiting the sh ort-period subdwa rf B of the eclipsing binary HW Virginis was first p resented in the works by 



planet system orbiting tne snort-period subdwart a ot tne eclipsing binary 11 W Virginis was hrst presented m tne works by 
I Lee et all (|2009] h and recentlv lBeuermann et all (|2010l ). |Potter et alj (|2011 ). and lDovle et al. 1 ()201 lh announced the existence 



of two circumbinary planets in possible mean motion resonances around NN Serpentis, two giant planets orbiting the eclipsing 
binary UZ Fornacis, and a single circumbina ry pl anet around the sta rs of the binary system Kepler 16, respectively. 

The studies bv lBeuermann et al. | l|20ld ) and iPotter et all |201ll ) inferred the presence of additional massive objects by 
explaining the observed timing anomalies with the light-travel time (LTT, hereafter) effect. In the ideal case, the stellar com- 
ponents of an eclipsing binary system are orbiting their common centre of mass with a constant period. Timing irregularities 
can occur if the binary system is accompanied by an additional massive object. In this case, the binary centre of mass is 
orbiting the system's centre of mass. At some times the binary will be closer to the observer while at other times, it will be 
farther away, giving rise to the LTT effect on measured eclipse egress times. 

In a recent study, the measurements of eclipse egress times of the eclipsi ng polar binary H U Aquarii (HU Aqr, hereafter) 
was used to infer the presence of a circumbinary planet around this system (|Qian et al.ll201ll h 

These authors modeled the complete timing data set by adding the LTT effects from two circumbinary planets and found 
that the pericentre of the outer planetary companion is inside the orbit of the inner planet. This orbita l architecture i mplies 
a crossing orbit configuration which points to strong mutual perturbations between the two planets. I Horner et alj (|201ll ) 
subsequently carried out a detailed stability analysis of these bodies and found that almost all their initial conditions lead to 
unstable orbits on short-time scales. These authors con cluded that the HU Aqr planetary system has most likely a different 
orbital architecture than proposed bv lQian et al.l (|201ll h 

Three possibilities exist to explain the instability of the proposed planetary syste m. Either i) the LT T parameter space 
was not explored thoroughly while orbital stability constrains were imposed [omitted in lQian et all (|201ll )]. or, ii) the applied 
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LTT model was not complete and was unable to explain the timing data set properly, or iii) the data set is not yet large 
enough to draw reliable conclusions on the existence of additional low-mass companions. 

In this work we aim to find a more plausible orbital architecture that best describes the timing data while also being con- 
form with orbital stabi lity requirements. We regard the stability co ndition as an observable that places additional constraints 
on the fitting process (|Gozdziewski. Konacki fc Macieiew ski 2005). We carried out a large-scale Monte Carlo, least-squared 
parameter survey by fitting various LTT models to the complete data set. In parallel to the fitting process, we performed 
a stability analysis of various LTT model orbits which best describe the data set. By parameterising planet's Hill radii, we 
imposed orbital stability constrain ts to our paramet er survey in order to explore orbital architectures that result in stable 
orbits. We extended the models bv lQian et al.l |2011) to systems that also allow the inner planet to attain an eccentric orbit. 
Our results suggest that two-planet LTT models in general tend to produce orbits with higher eccentricities which are likely 
unstable. 

The outline of this paper is as follows. In section 2 we present details of the adopted LTT model and provide a brief 
description of the least-squared minimisation algorithm. In section 3 we introduce our orbital stability constrains and in 
section 4 we describe results of numerical experiments. Finally, we give a summary and discussion in section 5. 



2 ANALYTIC LTT MODEL, LEAST-SQUARED FITTING WITH ORBITAL STABILITY 
CONSTRAINTS 

At the basis of our analysi s is th e combined ecli p se eg ress t iming data s et (a long with la error bars) as compiled from 
the works of ISchwope et all (|200ll ). TSchwarz et al.l ((2009), and lQian et al.l (|201 lh ■ A total of 113 timing measurements were 
obtained for HU Aqr which span the time between April 1993 (BJD 2449102.9) and May 2010 (BJD 2455335.3), equivalent to 
a period of approxim ately 17 years. All timing measurements have been stated in barycentric julian date (BJD) as described 
m iQian et all |201ll ). and thus have been transformed to the solar system barycentre. In the absence of a mechanism that 
causes timing variations, the linear ephemeris is given by 

BJD = T + P E, (1) {?} 

where E denotes the ephemeris cycle number, To is the reference epoch, and Po measures the eclipsing period of HU Aqr. 
The long observing baseline renders the eclipsing period to be known with high precision. In this work, we chose to place the 
reference epoch close to the middle of the observing baseline to avoid parameter correlations between To and Po during the 
fitting process. In the following, we outline our least-squared fitting technique, and develop the LTT model that describes the 
timing measurements by two independent planets, while imposing orbital stability constrains. 



2.1 Analytic LTT Model 

We implemented the LTT model as described in llrwinl 1 19521 ) in IDL and the resulting code is available upon request. In this 
model, the stars of the binary are assumed to represent one single object with a total mass equal to the sum of the masses of 
the two stars. The assumed single object is then placed at the original binary barycentre and in an orbit around the system's 
total centre of mass. In the following, we denote the latter as the LTT orbit. The reference system of the model originates 
from the centre of mass of all involved bodies with the z-axis pointing towards the solar system's barycentre. The quantity 
O — C for HU Aqr was then calculated using 

0-C = BJD oba -BJD = ^ Tk (a b sin I k , e b ,u b ,T b ,P b ), (2) {?} 

fc=l,2 

where BJD obs denotes the timing measurement series and, Ti,T2 are given by 



Tk 



2 



a b sinl k l — e b 



■ sin(/b + Ub) + e b sin tot 



l + e b cos f b 

In Eq.|3] the quantity c is the speed of light and f b is the true longitude of the binary's centre of mass in the LTT orbit. The 
LTT orbital parameters are the projected semi-major axis of the binary barycentre about the total centre of mass a b sin Ik, 
the binary eccentricity e b , argument of pericentre uj b , time of pericentre passage T b , and the period of the LTT orbit P b . The 
argument of pericentre is an inertial angle between the line of apse and a line that is perpendicular to the line connecting the 
solar system barycentre and the barycentre of the system (see fig. 1 in Irwin (1952)). We added a subscript to the inclination 
as it appears customary and is frequently used in the literature. However, both objects are orbiting in the same LTT plane. 
Therefore, they (the merged binary components and the unseen companion) share the same inclination. In general, the line- 
of-sight inclination of the LTT orbit with respect to the plane of the sky (Ik) is not known from the timing measurements. 
This is true except for cases where the binary companion reveals itself by eclipsing the binary components. In this case, the 
orbits are co-planar (e.g., the circumbinary planetary system of Kepler 16) and only minimum separation and minimum mass 
of the companion are obtained. 

To infer the instantaneous angular position of the bin ary's centre of mass requires the computation of the eccentric anomaly 
E' from the mean anomaly M (|Murrav fc DermottfcOOlh . Adopting the ephemeris cycle number (E) as the independent (time) 
variable, the mean anomaly is given by M — n(BJD — T) where n = 2ir/P is the mean-motion of the combined binary in its 
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LTT orbit. The eccentric anomaly is obtained from the orbital eccentricity and mean anomaly by solving Kepler's equation. 
For the purpose of this computation, we implemented the Newton-Raphson algorithm to iterate towards a solution for E' 
with an initial starting guess of E' — M. The true anomaly is then computed using 

We tested t he correct im plementation of the above model by reproducing the light-time curves for different values of ujb as 
presented in llrwinl l| 19591 ). 



2.2 Least-squared Fitting with MPFIT 

For the non-linear, least-squared fitting pro cess, we used the Leven berg-Marquardt (LM, hereafter) minimisation algorithm as 
implemented in the IDL routine MPFI tFI dMarkwardt et alj|2009h. This algorithm has found widespread applications within 



astronomy including timing-modeling (jBarlow. Dunlap fc Clemens! l201l] ). In this routine, partial derivatives for computing 
the gradient field in a parameter space are calculated from numerical differentiations of model functions. The routine also 
allows for imposing upper and lower boundary conditions for each parameter of the model. Parameters can be held fixed or 
flagged as variable as well. In short, MPFIT attempts to minimise the sum of weighted residuals squared (x 2 ) between the 
model and the observed timing measurements. The reduced x 2 is introduced as \r an d is given by x 2 divided by the number 
of degrees of freedom. The algorithm attempts to iterate towards a fit until the difference between two consecutive Xr are 
smaller than some user-defined threshold. 



2.3 The LTT Mass- Function of the Unseen Companion 

In essence, applying the LTT model to an anomalous timing data set provides information on the binaries barycentre motion 
and in principle no information is obtained about the unseen (too dim) companion. However, in the barycentric motion, the 
orbital period and eccentricity (in addition to the inclination of the LTT orbit) are shared between the close binary system 
and the unseen companion. From Kepler's law of motion and similar to the case of a single-line spectroscopic binary, the 
minimum mass of the third companion within the context of LTT can be expressed as 

47T 2 2 

m k sml k = -Qp5 Kca ■ ( 5 ){ ? } 

In this equation, vrik (where k = 1, 2) is the mass of the (probable) unseen companion, G is the gravitational constant, 
K = ai,sin/fc/c is the semi-amplitude of the observed light-travel time curve, and a is the minimum semi-major axis of the 
unseen companion relative to the binary barycentre. Without additional information about the unseen companion it would 
not be possible to determine its semi-major axis. However, from LTT considerations, one can write the mass-function of the 
companion (s) as 



47r 2 (a i) sin7fc) 3 _ Aty 2 (Kc) 3 _ ml sin 3 I, 



f(m k )- — GP 2 - {mb + mk y ( 6 H ? > 

This mass-function can be obtained by fitting for K and P considering mt as the (combined) mass of the two stars. We note 
that the equation above is transcendental in m k . Assuming that mt, is known, the minimum mass, m k sin I k , of the unsee n 
companion can be determined via iteration. The mass of the binary stars (m,b) have been reported bv lSchwope et all (|201lh . 
The primary white dwarf component was found to be 0.80 (±0.04) Mq while the secondary has a mass of 0.18 (±0.06) Mq 
at an inclination (determined from eclipses) of 87jl2'| degrees. Here we use rrib = 0.98 Mq. 

It is possible to show that the above transcendental equation has only one unique solution. In this work, we use the 
Newton-Raphson algorithm with /(rrifc) as the initial start guess. Finally, it is important to note the the LTT-determined 
semi-major axis and mass of the unseen circumbinary component are only calculated as projected quantities. 



3 LEAST-SQUARED FIT WITH ORBITAL STABILITY CONSTRAINTS 

In a similar spirit as outlined in iGozdziewski. Konacki fc Macieiewskil l|2003l . [20051 ). the dynamical stability of a multi-body 
system is treated as an "observable" . By "observable" we mean that a promising best-fitting model to the timing data should 
result in a stable orbit configuration as a necessary condition to substantiate its physical existence. 

To ensure the stability of the system, after each successful convergence towards a solution from the LM minimisation 
procedure, we tested the resulting fitted parameters for the following stability conditions: 

I. constraint: a± < a2 

II. constraint: 52 + nR^ > Qi + nR^ , with n = 1, 2, . . . 



1 http://purl.com/net/mpfit 
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where 02 = a2(l — £2) and Q\ = ai(l + ei) are the pericentre and apocentre for the outer and inner planet, respectively. The 
modified Hill radius for each planet is given by 

fl « = a .fi) 1/3 (7) 

and a multiple n of this quantity is used to parametrise the stability condition. It is important to note that is the semi- 
major axis of the two proposed planetary companions relative to the binary's centre of mass (still trea ted as a single objec t 
with mass mi,). At this point we note that we use a different definition of the Hill radius than used in iHorner et all l|201ll ). 
Our definition (although somewhat similar) results in a larger (around 40%) radius imposing a more stringent initial spacing 
condition between the two planets. 

From the two above-mentioned stability constraints, the first condition simply demands a hierarchical system with the 
size of the inner planet's orbit always being smaller than that of the outer body. For eccentric orbits, we use the "spacing 
parameter" n to quantify a stability condition. The second condition essentially demands that the outer planet's pericentre 
distance to be always larger than the inner p l anet's apocentre. This condition is clearly violated for the case of the planetary 
orbits in HU Aqr as proposed bv lQian et al.l l|201ll ). In an attempt to avoid close encounters, the condition is augmented by 
demanding that the spacing between the planets to be a multiple of their respective Hill's radii. For small values of n, mutual 
gravitational interactions are most likely to disrupt the system. Intuitively, for large values of n, the two orbits are more apart 
from each other. Hence imposing a larger spacing condition would constrain the parameter space considerably while possibly 
allowing for stable orbital architectures. To obtain stable orbits, after each iteration of the least-squared fitting procedure, 
we determine the minimum mass of the two companions. Combined with the orbital period, we then calculate the minimum 
semi-major axis asin/1,2 of the two planets from Kepler's third law. Finally, we evaluate the two stated stability conditions 
for certain orbits obtained from the fitting process. If a given set of fitted parameters passes the stability conditions, the 
resulting orbits are numerically integrated to directly infer their long-term stability. In this work we have used n = I as will 
be justified later. 



4 METHODS, NUMERICAL EXPERIMENTS, AND RESULTS 
4.1 Lomb-Scargle Frequency Analysis 

In this section, we describe our numerical experiments, the applied methods, and final outcome. As a first step, we applied a 
periodogram analysis to our dat a set using the P ERI0D04 algorithm. This analysis enables us to extract the frequency content 
in an unevenly spaced data set (|Lenz fe Br cgcr 2005) and also allows for an estimation of any periodicity that may exist in 
the data which can be used to provide an initial guess for the least-squared fitting procedure. At this stage, we are not aiming 
at accuracy. From the complete data set, we constructed the O — C diagram of HU Aqr at the measured times of eclipse 
egress, and used the period analysis algorithm. The resulting Lomb-Scargle powerspectrums are shown in Fig. [T] Two periods 
can clearly be seen in the data. A long-period variation of approximately 15.6 years, and a shorter period of 6.5 years. These 
values reflect the possible periods of two independent LTT orbits. Given the total observation time span of about 17 years, 
the shorter period is well-determined. This period has less power which implies a smaller orbit relative to the long period 
one. Later we will demonstrate that the timing data cannot be described by a single sinusoidal period alone. Therefore, both 
periodic signals must to some degree be real. 



4.2 Numerical Monte Carlo Experiments and Best-fitting Model 

In principle, searching for and finding a global minimum requires the probing of a large number of initial guesses spread across 
the Xr parameter space. Despite the vastness of this space, it is possible to limit the search by extracting additional information 
about some of the parameters of the system from the timing data set. For instance, the preceding period analysis can be used 
to obtain reasonable starting guesses for the fitting algorithm. Also, the model parameter at, sin 71,2 can be estimated from the 
largest amplitude in the data set. We find K ~ 3.47 x 10 -4 days as a maximum which translates into a projected semi-major 
axis of at, sin/2 = 0.06 au for the orbit of the outer planet. We carried out a "brute-force" search for a best-fitting model using 
a Monte Carlo approach. Initial guesses for the parameters of the model were then drawn from a random uniform and/or 
normal (Gaussian) distribution. The projected semi-major axis, argument of pericentre, and orbital eccentricities were chosen 
randomly from a uniform distribution with at, sin/1.2 G [0,0.060] au, uit G [— 7r,7r], and et, G [0,0.70]. We assigned an upper 
limit on the orbital eccentricities to explore different possibilities of orbital architectures while evaluating their goodness-of-fit 
to the data. The times of pericentre passages were also chosen at random from a uniform distribution with initial guesses 
bounded by the earliest and latest epochs of observation. Since we have a good idea of the initial reference epoch, eclipsing 
period, and the planet's orbital period, we decided to choose these values from a normal distribution. The mean and standard 
deviation for the reference epoch T was BJD 2452174.36592 and 0.0001 days, and for the orbital period P was 0.08682040612 
days and 10 -10 days, respectively. The mean of the LTT period Po was obtained from the Lomb-Scargle frequency analysis 
with standard deviation of two years for each planetary companion. In order to avoid parameter correlations between the 
refere nce epoch and the eclipsing period, we chose the reference epoch (E = 0) to be in the centre of the data set. lQian et al.l 
(2011) chose to place the reference epoch in their analysis at the first timing record. 
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We modeled the complete data set (as obtained from the literature) using a two-Keplerian LTT model. In their analysis 
and by construction. [Oian et al.l |201l[ ) forced the inner companion to be in a circular orbit. While this requirement reduces 
the number of fre e parameters, it does not account for the possibility of eccentricity excitation of the inner planet by an 
outside perturber (Malmberg fc Daviesll20oH '). In our study, we allowed the inner planet to also attain some eccentricity. It is 
however, important to point out that in the two-Kepler LTT model, no mutual gravitational perturbations are included. In 
our first approach to find a best-fitting model, we generated 111,844 initial guesses. Each LM guess was allowed a maximum 
of 500 iterations prior to termination. We carried out the Monte Carlo simulations on a supercomputer lasting for about 14 
days of non-stop computing. In detail, we used the two-sided numerical derivative option available in MPFIT algorithm. We 
also chose to assign individual maximum step-sizes to each model parameter. For example, we allowed the eccentricity to only 
change by at most 1% per iteration. Maximum change in an angular quantity was 5°. Similar reasonable values were assigned 
to other parameters. We believe that our initial choices of the value of the parameters renders the search-space of Xr to have 
a fine enough grid for each iteration of the fitting algorithm. We also believe that with such a large number of initial guesses 
the full Xr space will be sufficiently explored. During the fitting process, we held all twelve parameters of the system free, and 
decreased the numerical accuracy parameters in MPFIT from their default values. We also decreased the parameter step-size 
with which the numerical derivatives were evaluated. 

We then carried out a Monte Carlo experiment based on a total of 111844 randomly distributed initial guesses. For each 
converged initial guess with Xr < 10.0, we recorded the resulting Xr, the number of iterations, the initial guess parameters, 
and the final converged model parameters as returned by the LM minimisation algorithm. In order to avoid any bias effect, 
we removed converged iterations that reached a lower or upper boundary value. That means, once a parameter reached any of 
the boundaries within its defined range (see above), its value was held constant and the parameter was no longer changed by 
the minimisation algorithm for carrying out the remaining of the iterations. Although we have no strict explanation for this 
"stickiness" effect, we judge that it has the potential to skew the final distribution of model parameters. The idea is that from 
a mathematical point of view, the minimisation algorithm attempts to cross the lower or upper boundary to obtain a better fit 
based on the steepest descent in parameter space. Hence it remains constant at these locations but with no physical meaning. 
As an example, one can consider the orbital eccentricity. Since in the systems of our interest (circumbinary planetary systems) 
eccentricities smaller than zero do not correspond to acceptable orbits, iterations that were stuck at the lower boundary of 
e = were removed. This allows us to avoid introducing a bias towards a preferred direction in parameter space that has no 
reasonable physical significance. In addition, MPFIT algorithm does not return any formal uncertainties for fitted parameters 
that reach a lower or upper boundary. After this cleaning process, we ended up with a data file containing a total of 48485 
initial guesses. 

Our best fit resulted in a Xr — 1-43 which had an occurrence frequency of fourteen times in total. We show the cor- 
responding O — C diagram of this case in Fig. [2] Fitted model and derived parameters are shown in Table [1] with formal 
uncertainties as obtained from the solutions covariance matrix. We point out the close agreement between the two fitted 
periods (6.5 and 15.5 years, respectively) and the periods obtained from the Lomb-Scargle frequency analysis shown in Fig. [I] 
We visually inspected the remaining thirteen models and did not find any differences. All final fitte d parameters agre ed with 
one another. However, we noticed a qualitative difference when comparing our best fit with Fig. 2 in lQian et all (|201ll ). These 
authors show a better description of the newly added timing data but they do not provide a quantitative Xr value for their 
best fit. When comparing our results, we noticed that the inner compan ion in our model can attain eccentric orbits, which 
might explain the difference between the fits presented in this work and in lQian et al.l ({2011) . In addition, in Fig. [21 it appears 
that a single LTT (sinusoidal) curve will provide a poor description of the timing data. Neither the individual short (inner 
planet) nor the long-period (outer planet) LTT model will be able to provide a satisfactory fit by itself. Only by combining 
the two LTT models one can obtain a satisfactory fit. 

4.3 Orbital Integration of Best-fitting Model 

Our final orbital parameters for the two planetary companions are relatively close to ones reported by lOian et all (|201lh 
except for a slight increase (decrease) in the inner (outer) companion's orbital eccentricity. The outer companion's mass and 
semi-major axis are also larger in ou r work. We examined the orbital stability of o ur best-fitting model (Table [TJ using the 
orbital integration package MERCUR^ ( Chambers fc M igliorini 1997; Cha mbers! 1999T ). We integrated the three-body equations 
of motion using the variable-step Bulirsh-Stoer (BS2) algorithm. Initial step-size was chosen to be 1% of a day and the accuracy 
parameter was set to 10 -16 . To be consistent with the fitted LTT orbits, we considered the binary pair as one single object 
with a mass equal to the sum of the masses of the two stars. We carried out integrations for two scenarios: i) the nominal 
parameters for the best-fitting model with (xr = 1-43), and ii) an optimistic scenario where we considered the lower limit 
in masses and eccentricity and upper limit in the semi-major axis for both companions. We integrated both systems for 10 4 
years. The results are shown in Fig. [3] This short integration time is sufficient to demonstrate an unstable configuration for 
the first scenario. The upper panel of the figure shows the orbits based on the nominal best-fitting parameters and clearly 
shows orbital instability of at least one of the planet companion that resulted in its eventual escape. 

The optimistic scenario shown in the lower panel of Fig. 3 resulted in a relatively stable orbital architecture when compared 
with the nominal scenario. The corresponding Kepler elements of the orbit in this case and their time evolutions are shown 



2 www.arm.ac.uk/~jec 
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in Fig. [4] From this figure we note that both our initial stability conditions are satisfied: i) the inner companion has a smaller 
semi-major axis than the outer planet and, ii) the quantity q2 + R^p ~ 5.9 au stay always larger than Qi + Rff ~ 4.2 au. 
Also in this scenario, the two companions are well separated from one another. The latter justifies our choice of n — 1 for the 
stability spacing parameter as introduced earlier. As shown in the top panel of Fig. 3, the parameters of our nominal best- 
fitting model resulted in an unstable orbital configuration. Combined with the results of the optimistic scenario experiment, 
this suggests that stable systems are to be found among less eccentric orbits for the inner and/or outer companion. 



4.4 Xr Frequency Distribution 

In order to find model parameters that result in a stable orbital configuration, we closely examined the frequency distribution 
of the Xr values. Prior to our experiment, we had expected a (l/Xr)-like histogram distribution. In this distribution the 
majority of final Xr values will belong to the first bin with smallest Xr- ln Fig. [5]we show the histogram with a bin size of 0.10 
that was obtained from our actual results. As is evident from this figure, increasing the bin size would result in the expected 
1/Xr distribution. However, we noticed that unlike what is expected from this type of distribution, a non-negligible fraction 
(4%) of initial guesses produce models with final fits in the interval 1.43 ^ Xr < 1-53. The histogram's peak value is associated 
with the interval 1.53 ^ Xr < 1-63 amounting to 10% of the total number of guesses. Technically, our best-fitting model with 
X r = 1.43 represents the global minimum in the Xr _ P arame t er space, and hence gives the most precise description of the data. 
Though, the majority of initial guesses have a higher frequency at somewhat larger Xr values. 

We also examined the distribution of the final fitted eccentricities (e) and minimum semi-major axes (a^sin/i^) for 
the two companions for all bins shown in the histogram of Fig. [S] In Fig. \§\ we show the resulted distribution for the final 
eccentricities and in Fig. [7] for the final semi-major axes. For the first bin interval (Fig. [6^,) starting with Xr = 1-43, our 
simulations indicate that the final fitted eccentricity for the outer companion is in general high with a minimum at around 
0.20 and a maximum at 0.50. The final eccentricity for the inner planet ranges from sa 0.0 to 0.30. The corresponding final 
fitted semi-major axes appear to be less scattered with final fitted values forming an arc-like feature as shown in Fig. [7] In 
particular, we note a lower limit on the projected semi-major axis of the outer planet equal to 0.0312 au. The inner companion's 
semi-major axis appears to be constrained between 0.008 and 0.02 au. The cross hair in the upper left panel in Figs. [5] and [7] 
represents the location of the average of all fourteen best-fitting models (with Xr = 1-43) final minimum semi-major axes and 
eccentricities. Based on our previous orbital stability analysis of the nominal best-fitting orbits, we estimate that the final 
fitted parameters obtained from the first Xr bin, will result in unstable orbits because of large orbital eccentricities that cause 
close approaches between the two planets. 



4.5 New System Parameters and Orbital Stability Constraints 

For the first three bin intervals in Fig. [51 the final eccentricities scatter over a larger range. The majority of the initial parameter 
guesses (4842 in total) ended up with Xr between 1.53 and 1.63. Our results point to a qualitative emergence of an apparent 
bimodal distribution for the bin interval 1.83 ^ Xr < 1-93. A density enhancement is observed for final eccentricities centered 
at (0.07,0.3) and a second one at (0.12,0.01). Although these final fitting parameters correspond to a poorer goodness-of-fit 
quality, they allow near-circular orbits for the outer companion and hence would satisfy our stability constraints as outlined 
earlier. 

An inspection of Fig. 7 indicates that the bin interval 1.93 ^ Xr < 2.03 also seems to result in circular orbits for 
both companions, showing an enhancement of final eccentricities around (0,0). We would like to stress that these findings 
are of qualitative judgment from visual inspections. We tested the orbital stability of one of the final parameters with 
1.83 ^ Xr < 1-93. The final eccentricity of the two planetary orbits were close to circular. We have shown the locations of 
these orbits with a haircross in Figs. |6^ and [7^ for their corresponding semi-major axes. The associated fitted model had 
Xr = 1-89 and is shown in Fig. [8] Table [2] shows the corresponding numerical values of the system's parameters. Comparing 
these results with those shown in Fig. [5] indicates that the qualitative differences are minute but with dramatic consequences 
for the underlying stability of the proposed circumbinary planetary system. 

We also tested the orbital stability of the model with Xr = 1-89. The result of a 1000-year integration is shown in Fig. [9] 
We show the pericentre, semi-major axis, and apocentre distances in the top panel and the eccentricity in the bottom panel. 
The fitted parameters in this system satisfy the two stability conditions outlined before. While the semi-major axis exhibits 
short-term variations, the variations in the two eccentricities are long-term and secular. This is possibly due to exchange of 
orbital angular momentum between the two planetary companions. We also carried out the integration over a 10 6 year time 
interval and found the orbits to be confined and stable with no signs of global instabilities. Particularly interesting is that the 
Xr = 1-89 fit parameters do not show an orbit-crossing architecture. To substantiate our results, we examined several other 
orbital parameters with goodness-of-fits in the neighbourhood of Xr = 1-89. All of the final parameters resulted in stable 
orbits over a 10 6 year integration period. Carrying out a detailed dynamical analysis of the orbits shown in Fig. [9] is beyond 
the scope of this work. 
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5 SUMMARY AND DISCUSSION 

In a recent study iQian et al l l|201ll ) announced the discovery of a two-planet system around the eclips ing polar binary HU 
Aquarii. In a subsequent paper, the stability of their best-fitting orbital parameters was examined by Irlorner et all (l201lf ) 
who showed that the system is highly unstable. The most likely reason for this instability is the high eccentricity of the outer 
companion. In an attempt to find a system with stable orbits, we calculated a large number of two-planet LTT models using 
a Monte Carlo approach. This approach provides statistics on the distribution of the final goodness-of-fit Xr parameter in 
combination with constrains on requiring initial orbital stability. 

In this work, we analyzed the complete available timing data of HU Aquarii. We made the following assumptions in 
our analysis. First we allowed the inner planet to attain an e c centri c orbit in response to possible eccentricity excitation 
by an outside perturber. Compare to the study by IQian et al.l (|201ll ). this assumption resulted in introducing three extra 
free parameters. We also did not include a secular term which would account for timi ng variations due to mass-/angular 
momentum transfer and/or magnetic activity of the binary component. IQian et alj (|201lh argue that the secular term in their 
work is most likely caused by the presence of a third circumbinary massive object. However, the existence of this body in 
their data is far from obvious as the observational time span is much smaller than the orbital period of the proposed third 
circumbinary (planetary) companion. The confirmation of the third planet would require a few decades of future photometric 
follow-up monitoring. Finally, we did not consider mutual gravitational interactions between the two planets or gravitational 
perturbations on the binary stars. 

Our best fit model had a goodness-of-fit parameter of Xr = 1.43 occurring fourteen times out of 111.844 initial guesses. 
From studying the histogram distribution of final Xr parameters, we were faced with a dichotomy. On the one hand, our 
best-fitting model of Xr = 1-43 had a low occurrence frequency (Fig. [5} and resulted in high eccentricities (Fig. [6ji). On 
the other hand, we encountered a significant higher occurrence frequency of larger Xr values. At the moment we have no 
explanation for this trend. One possibility is the inability of the LM minimisation algorithm to escape a local minimum in the 
non-linear parameter space. To test this possibility, one can use a least-squared minimisation procedure based on a genetic 
algorithm (GA) or a Bayesian Markov Chain Monte Carlo method. 

The results of our study suggest that we should be able to trust the nominal Xr = 1-43 best-fitting parameters as this fit 
represents, most likely, the global minimum in Xr parameter-space, and provides the best possible description of the observed 
timing data set. However, as we have shown, the resulting best-fitting orbit with Xr = 1-43 was unstable. Models in the 
neighbourhood of this system were also unstable with a collision or ejection taking place on very small time scales. This 
instability can be attributed to high values of the final fitted eccentricities for one or both planetary companions. It is worth 
mentioning that models with goodness-of-fit close to the nominal best-fitting model did not allow low-eccentricity /circular 
orbits which would have rendered them most likely to be stable. We would like to mention that other author s have encountered 
similar (in) stability pro blem for the two planets in the system s of NN Serpent is and UZ Fornacis as well (|Potter et alj|20ill : 
iBeuermann et al"1l2010l ). The two-planet system HW Virginis jLee et alj|2009l ) also appears to be unstable (B. Funk, private 
communication). In particular, it is interesting to note that the best-fittin g models of the proposed UZ Fornax circumbinary 
planetary system also attained high final eccentricities |Potter et al.ll201lT ). 

Examining the Xr frequency distribution in detail resulted in the emergence of a second population of models with fitted 
parameters that fulfill our orbital stability constraints. Orbital stability was examined for models with Xr = 1-89. Although 
this population had a larger Xr goodness-of-fit than the nominal case, we were able to find stable orbits due to significantly 
lower orbital eccentricities. Combining/pairing a least-squared minimisation algorithm with a stabilit y analysis enabled us to 
determine a new set of orbital elements for the two planets which are different from those given bv lQian et al.l (|201ll ). The 
latter is another implication of our study which points to the existence of a variety of stable configurations when considering 
models with increasingly larger Xr- 

Assuming that no other physical effects are capable of explaining the observed timing variations in the light curve of 
HU Aquarii, we are left with the two-planet model as the most probable cause for the observed timing variation. Our study 
indicates that single one-planet LTT sinusoidal timing variations are inadequate to describe the observational data. However, 
it is necessary to mention that timing anomalies could also be caused by direct perturbations of a circumbinary planet on the 
binary orbit, though these perturbations would introduce a high-frequency (short-period) component to the eclipse timing 
variations due to the short orbital period of the binary system. A second possibility worthy to consider in future work would be 
to include mutual perturbations between the two proposed circumbinary planets. Explaining eclipse timing variations caused 
by the latter effect has not been explored for circumbinary planetary systems. 

Finally, a similar study bv lWittenmver et alj (|201ll ) also addresses the instability of the two proposed companions around 
HU Aqr from revised LTT models. In their work the authors point out the need to obtain a better understanding of the mutual 
interaction between the two binary components. Binary orbital period modulation caused by magnetic interaction and/or 
mass-transport between the two binary components could also possibly explain the observed timing anomaly for HU Aqr. 
However, a combination of the above effects might also provide a satisfactory description of the observed timing data. At this 
point, we recommend future discoveries of multi-planet circumbinary systems also include a preliminary orbital stability study 
as a necessary condition to increase the likelihood of the existence of such systems. At the moment it seems that best-fitting 
single LTT models in superposition (n +T2) arc inadequate to result in stable multi-planetary circumbinary systems. Though 
such a model provide a good description to the observational timing data. Furthermore, we also recommend future studies 
to obtain and publish timing data of the secondary eclipse (if available). Any timing variations introduced from gravitational 
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interactions and/or LTT should also provide an adequate description of timing measurements obtained from the secondary 
eclipse. 

Further constraining of our model requires obtaining additional timing data of the HU Aqr system. Using the results of 
our study involving stable planetary orbits, we predict that future mid-egress timing measurements would result in O — C 
timing differences of -10.0 to 0.0 seconds (predicted at epoch E = 40000 in Fig. [8]). We emphasise that high-accuracy timing 
measurements are crucial to unveil the true nature of the observed timing anomalies as well as the development of improved 
models describing eclipse (egress) timing variations of a binary star system. 
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Figure 1. Lomb-Scargle pcriodogram of the O — C eclipse egress timing measurements using the complete data set from the literature. 
Two distinct frequencies are found with /i = 0.000176 cycles/day (15.6 years) and fa = 0.000425 cycles/day (6.5 years). The amplitude 
of l//i is larger because it corresponds to the period of the wide orbit companion. 
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Figure 2. Best-fitting two-Kepler LTT model with reduced Xr = 1-43. The two single LTT sinusoids are superimposed. The inner 
companion is on a smaller orbit with a smaller am plitude. The smalle st and largest timing errors are 0.10 and 9.99 seconds, respectively. 
The group of data at cycle number 34000 are from lQian et al. I fcoiily 



Table 1. Best-fitting parameters for the LTT orbits of HU Aqr. Subscript b refers to orbital quantities for the compined binary. Subscript 
1, 2 refers to the circumbinary companions. Note that the eccentricity, inclination and orbital period are shared quantities as outlined in 
the text. The last two entries denote the minimum mass and minimum semi-major axis of the two proposed companions. 



Parameter 


two-Kepler LTT 


Unit 


ri T2 


red. Xr 


1.43 





To 


2,452,174.365694(2) 


BJD 


Po 


0.0868204066(14) 


days 


a b sinii >2 


0.0179(5) 


0.0353(5) 


au 


e& (or ei,2) 


0.075(35) 


0.29(3) 






316(2) 


351(1) 


deg 


n 


2,453,462(16) 


2,450,050(7) 


BJD 


P b (orPi.a) 


2359(7) 


5646(8) 


days 


K 


1.03(0.03) x 10" 4 


2.04(0.02) x 10" 4 


days 


f(M 1>2 ) 


1.38(0.13) x 10~ 7 


1.84(0.09) x 10~ 7 


M Jup 


Mi >2 sin/i^ 


5.4(0.2) 


5.9(0.1) 


M Jup 


ai : 2 sin 71,2 


3.45 


6.18 


au 
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Figure 3. Orbital stability of the two companions. Upper panel: Unstable system corresponding to the nominal case of best-fitting 
(Xr = 1-43) model. Lower panel: Stable system for the optimistic parameters (see text). 
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Summing all bins results in 48485 guesses (see text). 



10.0. 



Table 2. Best-fitting parameters for the LTT orbits of HU Aqr corresponding to Fig. \E\ Subscript b refers to orbital quantities for 
the compined binary. Subscript 1, 2 refers to the circumbinary companions. Note that the eccentricity, inclination and orbital period 
are shared quantities as outlined in the text. The last two entries denote the minimum mass and minimum semi-major axis of the two 
proposed companions. 



Parameter 


two-Kepler LTT 


Unit 




Tl 


T2 




red. Xr 




1.89 




To 


2,452,174.3657082(5) 


BJD 


Pq 


0.0868204062(10) 


days 


a b sin7i i2 


0.0146(7) 


0.0321(5) 


au 


e& (or ei i2 ) 


0.11(3) 


0.04(3) 






292(2) 
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deg 


n 


2450992(15) 


2449837(6) 


BJD 
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2226(9) 


5155(8) 


days 


K 


8.43(0.04) x 10- 


5 1.85(0.01) x 10" 4 


days 


/(Mi, 2 ) 


8.38(0.11) x 10" 


8 1.66(0.07) x 10- 7 


M Jup 


Mi ;2 sin/i j2 


4.5(0.2) 


5.7(0.2) 


M Jup 


01,2 sin7i. 2 


3.32 


5.81 


au 
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Figure 6. Scatter plots of final fitted eccentricities for different x^-bin intervals following the histogram distribution shown in Fig. \E\ 
An apparent bimodal distribution is found for the bin 1.83 ^ Xr < 1-93. We note that all orbits have 0.0 < ei.2 < 0.70 since we removed 
those initial guesses that "stick" to the lower and upper boundary of final eccentricity during the minimisation/iteration process. 
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Figure 7. Same as Fig. [6] but this time for the final semi-major axes of the two proposed planetary companions. The semi-major axis 
is the fitted minimum semi-major axis sin/1,2. 



Least-squared modeling with orbital stability constraints on planets in HU Aquarii 17 




Figure 9. Result of a 1000-year orbital integration of the two proposed companions for the best-fitting initial condition with \r = 1-89. 
The integration was extended to 10 6 years and showed bounded/stable motion. Upper panel: Time evolution of periastron, semi-major 
axis and apastron distance. Bottom panel: Time evolution of eccentricity. 



